program uva_amp_gen
    use, intrinsic :: iso_fortran_env, only: real64
    implicit none
    ! 定义参数
    integer      :: nh ! 谐波次数
    real(real64) :: t0 ! 激励起始时间 (s)
    real(real64) :: omega! 基频角速度 (rad/s)
    real(real64) :: a0 ! 静态偏移量
    real(real64),allocatable :: an(:) ! 余弦项振幅
    real(real64),allocatable :: bn(:) ! 正弦项振幅
    real(real64) :: t_start ! 时间范围起始
    real(real64) :: t_end ! 时间范围结束
    real(real64) :: dt ! 时间步长
    integer :: n_points ! 数据点数量
    character(len=200)::result_file
    ! 声明变量
    real(real64), dimension(:),allocatable :: t, tau, u, v, a
    integer :: i, n,fn
    logical, dimension(:),allocatable :: mask                ! 用于标记t >= t0的区域

! -----------------------------
!        读取输入
! -----------------------------
    print *, 'You nees define this formula: '
    print *, 'A(t)='
    print *, '     (t_start <= t <= t_end)'
    print *, '     when t<t_0, then A_0'
    print *, '     when t>t_0, then A_0 + A_n*cos(n*\omega*(t-t_0)) + B_n*sin(n*\omega*(t-t_0)) [n=1~N]'
    print *, ''
    print *, 'Input harmonic number [N] ='
    read(*,*) nh
    allocate(an(nh),bn(nh))
    print *, 'Input start time [t_0] ='
    read(*,*) t0
    print *, 'Input base frequency angular velocity [\omega] ='
    read(*,*) omega
    print *, 'Input Static offset [A_0] ='
    read(*,*) a0
    print *, 'Input A_n values, use whitespace to delimeter, need ',nh,' floats'
    read(*,*) an
    print *, 'Input B_n values, use whitespace to delimeter, need ',nh,' floats'
    read(*,*) bn
    print *, 'Input t_start ='
    read(*,*) t_start
    print *, 'Input t_end ='
    read(*,*) t_end
    print *, 'Input Time step [dt] ='
    read(*,*) dt
    print *, 'Input Amplitude curve data filename ='
    read(*,'(a)') result_file
    
    n_points= int((t_end - t_start)/dt) + 1  
    allocate(t(n_points))
    allocate(tau(n_points))
    allocate(u(n_points))
    allocate(v(n_points))
    allocate(a(n_points))
    allocate(mask(n_points))
    ! 生成时间数组
    do i = 1, n_points
        t(i) = t_start + (i-1)*dt
        tau(i) = t(i) - t0
        mask(i) = (t(i) >= t0)
    end do

    ! 计算位移u(t)
    u = a0  ! 初始化为静态偏移a0
    do n = 1, nh
        do i = 1, n_points
            if (mask(i)) then
                u(i) = u(i) + an(n)*cos(real(n, real64)*omega*tau(i)) + &
                            bn(n)*sin(real(n, real64)*omega*tau(i))
            end if
        end do
    end do

    ! 计算速度v(t)（u的一阶导数）
    v = 0.0_real64
    do n = 1, nh
        do i = 1, n_points
            if (mask(i)) then
                ! 余弦项导数：d/dt [An*cos(nωτ)] = -An*nω*sin(nωτ)
                ! 正弦项导数：d/dt [Bn*sin(nωτ)] = Bn*nω*cos(nωτ)
                v(i) = v(i) - an(n)*real(n, real64)*omega*sin(real(n, real64)*omega*tau(i)) + &
                            bn(n)*real(n, real64)*omega*cos(real(n, real64)*omega*tau(i))
            end if
        end do
    end do

    ! 计算加速度a(t)（u的二阶导数）
    a = 0.0_real64
    do n = 1, nh
        do i = 1, n_points
            if (mask(i)) then
                ! 余弦项二阶导数：d²/dt² [An*cos(nωτ)] = -An*n²ω²*cos(nωτ)
                ! 正弦项二阶导数：d²/dt² [Bn*sin(nωτ)] = -Bn*n²ω²*sin(nωτ)
                a(i) = a(i) - an(n)*(real(n, real64)**2)*(omega**2)*cos(real(n, real64)*omega*tau(i)) - &
                            bn(n)*(real(n, real64)**2)*(omega**2)*sin(real(n, real64)*omega*tau(i))
            end if
        end do
    end do

    open(newunit=fn,file=trim(result_file),action='write',status='replace')
    do i=1, n_points
        write(fn,'(4F22.15)') t(i),u(i),v(i),a(i)
    end do
    close(fn)
end program uva_amp_gen